Chapter 23: Model basic

23. Model basics

23.1 Introduction

Trong model sẽ bao gồm 2 phần cần chú ý:

1. Family: là dạng của model muốn sử dụng, tuyến tính hay hàm số bậc hai, ba
2. Fitted model: để fit dạng family model muốn sử dụng gần đúng nhất với data

Quan trọng nhất là model “tốt nhất” (theo tuỳ chuẩn) không có nghĩa là model đúng với thực tế “true”

All models are wrong, but some are useful. - George Box

For such a model there is no need to ask the question “Is the model true?”. If “truth” is to be the “whole truth” the answer must be “No”. The only question of interest is “Is the model illuminating and useful?”.

Đầu tiên load các libary cần thiết

library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.4.1     ✔ purrr   1.0.1
## ✔ tibble  3.1.8     ✔ dplyr   1.1.0
## ✔ tidyr   1.3.0     ✔ stringr 1.5.0
## ✔ readr   2.1.3     ✔ forcats 1.0.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()

library(modelr)
options(na.action = na.warn)

23.2 Model cơ bản

Đầu tiên xem data set của sim1 trong model r, gồm 2 biến x và y ở dạng continuous

sim1

## # A tibble: 30 × 2
##        x     y
##    <int> <dbl>
##  1     1  4.20
##  2     1  7.51
##  3     1  2.13
##  4     2  8.99
##  5     2 10.2 
##  6     2 11.3 
##  7     3  7.36
##  8     3 10.5 
##  9     3 10.5 
## 10     4 12.4 
## # … with 20 more rows

ggplot(sim1, aes(x, y)) +
         geom_point()

Từ plot có thể thấy hai biến x và y có correlation mạnh và ở dạng tuyến tính, y = a_0 + a_1*x.

Thử dùng geom_abline để vẽ các đường thẳng với intercept và slope bất kỳ

models <- tibble(
  a1 = runif(250, -20, 40),
  a2 = runif(250, -5, 5)
)

ggplot(sim1, aes(x, y)) + 
  geom_abline(aes(intercept = a1, slope = a2), data = models, alpha = 1/4) +
  geom_point() 

=> Phải tìm giá trị a_0 và a_1 khớp nhất với data.

Sử dụng khoảng cách khác biệt giá trị y dự đoán và giá y trị thực tế để tìm model thích hợp nhất.

Ví dụ, tính giá trị y dự đoán nếu a_0 = 7 và a_1 =1.5

y = a_0 + a_1 *x

y = 7 + 1.5 *x

model1 <- function(a, data) {
  a[1] + data$x * a[2]
}
model1(c(7, 1.5), sim1) #y = 7 + 1.5 *x 

##  [1]  8.5  8.5  8.5 10.0 10.0 10.0 11.5 11.5 11.5 13.0 13.0 13.0 14.5 14.5 14.5
## [16] 16.0 16.0 16.0 17.5 17.5 17.5 19.0 19.0 19.0 20.5 20.5 20.5 22.0 22.0 22.0

viết function để tính độ lệch

measure_distance <- function(mod, data) {
  diff <- data$y - model1(mod, data)
  sqrt(mean(diff ^ 2))
}
measure_distance(c(7, 1.5), sim1)

## [1] 2.665212

Dùng purrr package để tính căn bậc hai của sai số bình phương cho 250 models

sim1_dist <- function(a1, a2) {
  measure_distance(c(a1, a2), sim1)
}

models <- models %>% 
  mutate(dist = purrr::map2_dbl(a1, a2, sim1_dist))
models

## # A tibble: 250 × 3
##         a1      a2  dist
##      <dbl>   <dbl> <dbl>
##  1  -6.62   3.94    5.86
##  2  22.2    0.186   9.60
##  3  11.9   -0.0713  7.62
##  4  -0.692  2.61    3.23
##  5  14.8   -1.29   12.6 
##  6  32.0    1.86   26.9 
##  7   6.70   3.95   14.2 
##  8 -19.9   -1.88   47.2 
##  9   2.43  -2.86   32.1 
## 10 -12.6   -0.293  30.6 
## # … with 240 more rows

Thử vẽ 10 model với độ lệch nhỏ nhất

ggplot(sim1, aes(x, y)) + 
  geom_point(size = 2, colour = "grey30") + 
  geom_abline(
    aes(intercept = a1, slope = a2, colour = -dist), 
    data = filter(models, rank(dist) <= 10)
  )

Ngoài ra còn có thể chuyển giá trị a_1 a_2 vào scatterplot để xem

ggplot(models, aes(a1, a2)) +
  geom_point(data = filter(models, rank(dist) <= 10), size = 4, colour = "red") +
  geom_point(aes(colour = -dist))

Bỏ scatterplot này vào lưới để chọn giá trị a1 và a2 một cách đồng đều hơn

grid <- expand.grid(
  a1 = seq(-5, 20, length = 25),
  a2 = seq(1, 3, length = 25)
  ) %>% 
  mutate(dist = purrr::map2_dbl(a1, a2, sim1_dist))

grid %>% 
  ggplot(aes(a1, a2)) +
  geom_point(data = filter(grid, rank(dist) <= 10), size = 4, colour = "red") +
  geom_point(aes(colour = -dist)) 

Đem 10 models theo scatter plot này vào lại data

ggplot(sim1, aes(x, y)) + 
  geom_point(size = 2, colour = "grey30") + 
  geom_abline(
    aes(intercept = a1, slope = a2, colour = -dist), 
    data = filter(grid, rank(dist) <= 1)
  )

Nếu làm cho grid càng nhỏ thì cuối cùng ta sẽ có model fit “best”.

Hoặc sử dụng phương pháp Newton-Raphson: chọn một điểm rồi tìm slope dốc nhất.

best <- optim(c(0, 0), measure_distance, data = sim1)
best$par

## [1] 4.222248 2.051204

ggplot(sim1, aes(x, y)) + 
  geom_point(size = 2, colour = "grey30") + 
  geom_abline(intercept = best$par[1], slope = best$par[2])

Phương pháp này ứng dụng được trên mọi dạng model, nhưng trong model linear, có function lm để fit model tốt nhất với data

sim1_mod <- lm(y ~ x, data = sim1)
coef(sim1_mod)

## (Intercept)           x 
##    4.220822    2.051533

23.2.1 Exercices

23.3 Vẽ model

23.3.1 Dự đoán

Tạo data dự đoán

grid <- sim1 %>%
  data_grid(x)
grid

## # A tibble: 10 × 1
##        x
##    <int>
##  1     1
##  2     2
##  3     3
##  4     4
##  5     5
##  6     6
##  7     7
##  8     8
##  9     9
## 10    10

Thêm giá trị y dự đoán vào

sim1

## # A tibble: 30 × 2
##        x     y
##    <int> <dbl>
##  1     1  4.20
##  2     1  7.51
##  3     1  2.13
##  4     2  8.99
##  5     2 10.2 
##  6     2 11.3 
##  7     3  7.36
##  8     3 10.5 
##  9     3 10.5 
## 10     4 12.4 
## # … with 20 more rows

#y = a_1 + a_2 * x
grid <- grid %>% 
  add_predictions(sim1_mod) 
grid

## # A tibble: 10 × 2
##        x  pred
##    <int> <dbl>
##  1     1  6.27
##  2     2  8.32
##  3     3 10.4 
##  4     4 12.4 
##  5     5 14.5 
##  6     6 16.5 
##  7     7 18.6 
##  8     8 20.6 
##  9     9 22.7 
## 10    10 24.7

Vẽ data dự đoán vào trong plot data quan sát

ggplot(sim1, aes(x)) +
  geom_point(aes(y = y)) +
  geom_line(aes(y = pred), data = grid, colour = "red", linewidth = 1)

Residuals

Function add_residuals() để thêm giá trị residual vào bảng prediction

sim1 <- sim1 %>% 
  add_residuals(sim1_mod)
sim1

## # A tibble: 30 × 3
##        x     y    resid
##    <int> <dbl>    <dbl>
##  1     1  4.20 -2.07   
##  2     1  7.51  1.24   
##  3     1  2.13 -4.15   
##  4     2  8.99  0.665  
##  5     2 10.2   1.92   
##  6     2 11.3   2.97   
##  7     3  7.36 -3.02   
##  8     3 10.5   0.130  
##  9     3 10.5   0.136  
## 10     4 12.4   0.00763
## # … with 20 more rows

Cách phương pháp vẽ graph cho residual

ggplot(sim1, aes(resid)) + 
  geom_freqpoly(binwidth = 0.5)

Hoặc ở dạng scatterplot

ggplot(sim1, aes(x, resid)) + 
  geom_ref_line(h = 0) +
  geom_point() 

Scatterplot của resid ở dạng random noise -> không có bias

plot(sim1_mod)

23.3.3 Exercises

Dùng model_matrix() để xem dạng model

df <- tribble(
  ~y, ~x1, ~x2,
  4, 2, 5,
  5, 1, 6
)
df

## # A tibble: 2 × 3
##       y    x1    x2
##   <dbl> <dbl> <dbl>
## 1     4     2     5
## 2     5     1     6

model_matrix(df, y ~ x1)

## # A tibble: 2 × 2
##   `(Intercept)`    x1
##           <dbl> <dbl>
## 1             1     2
## 2             1     1

dạng model là: y ~ x1 tương đương với y =

Thêm -1 để bỏ cột intercept mặc định

model_matrix(df, y ~ x1 - 1)

## # A tibble: 2 × 1
##      x1
##   <dbl>
## 1     2
## 2     1

y = 2 + x

model_matrix(df, y ~ x1 + x2)

## # A tibble: 2 × 3
##   `(Intercept)`    x1    x2
##           <dbl> <dbl> <dbl>
## 1             1     2     5
## 2             1     1     6

23.4 Formulas và model families

Model cơ bản y ~ x viết ở dưới dạng y = a_1 + a_2 * x

Để xem cách R làm model, dùng function model_matrix

df <- tribble(
  ~y, ~x1, ~x2,
  4, 2, 5,
  5, 1, 6
)
model_matrix(df, y ~ x1)

## # A tibble: 2 × 2
##   `(Intercept)`    x1
##           <dbl> <dbl>
## 1             1     2
## 2             1     1

R tự động thêm cột intercept vào model với cột đầu chỉ có số 1. Cách để bỏ cột số một trong model_matrix

model_matrix(df, y ~ x1 -1 )

## # A tibble: 2 × 1
##      x1
##   <dbl>
## 1     2
## 2     1

model_matrix(df, y ~ x1 + x2)

## # A tibble: 2 × 3
##   `(Intercept)`    x1    x2
##           <dbl> <dbl> <dbl>
## 1             1     2     5
## 2             1     1     6

23.4.1 Dạng categorical

Tạo data frame mẫu và model response phụ thuộc vào giới tính

df <- tribble(
  ~ sex, ~ response,
  "male", 1,
  "female", 2,
  "male", 1
)
model_matrix(df, response ~ sex)

## # A tibble: 3 × 2
##   `(Intercept)` sexmale
##           <dbl>   <dbl>
## 1             1       1
## 2             1       0
## 3             1       1

sim2 %>% head()

## # A tibble: 6 × 2
##   x         y
##   <chr> <dbl>
## 1 a     1.94 
## 2 a     1.18 
## 3 a     1.24 
## 4 a     2.62 
## 5 a     1.11 
## 6 a     0.866

ở dạng graph

ggplot(sim2) +
  geom_point(aes(x, y))

fit model vào

mod2 <- lm(y ~ x, data = sim2)

grid <- sim2 %>% 
  data_grid(x) %>% 
  add_predictions(mod2)
grid

## # A tibble: 4 × 2
##   x      pred
##   <chr> <dbl>
## 1 a      1.15
## 2 b      8.12
## 3 c      6.13
## 4 d      1.91

Prediction của y ở mỗi x chính là trung bình của y ở điểm đó

ggplot(sim2, aes(x)) + 
  geom_point(aes(y = y)) +
  geom_point(data = grid, aes(y = pred), colour = "red", size = 4)

Interaction (continuous and categorical)

Trong trường hợp sim 3 vừa có continuous and categorical

ggplot(sim3, aes(x1, y)) + 
  geom_point(aes(colour = x2))

2 models có thể fit cho data này

mod1 <- lm(y ~ x1 + x2, data = sim3)
mod2 <- lm(y ~ x1 * x2, data = sim3)

(y ~ x1 + x2): y = a_0 + a_1*x1 + a_2*x2

(y ~ x1 * x2): y = a_0 + a_1*x1 + a_2*x2 +a_12*x1*x2

Tạo bảng prediction

grid <- sim3 %>% 
  data_grid(x1, x2) %>% 
  gather_predictions(mod1, mod2)
grid

## # A tibble: 80 × 4
##    model    x1 x2     pred
##    <chr> <int> <fct> <dbl>
##  1 mod1      1 a      1.67
##  2 mod1      1 b      4.56
##  3 mod1      1 c      6.48
##  4 mod1      1 d      4.03
##  5 mod1      2 a      1.48
##  6 mod1      2 b      4.37
##  7 mod1      2 c      6.28
##  8 mod1      2 d      3.84
##  9 mod1      3 a      1.28
## 10 mod1      3 b      4.17
## # … with 70 more rows

Graph

ggplot(sim3, aes(x1, y, colour = x2)) + 
  geom_point() + 
  geom_line(data = grid, aes(y = pred)) + 
  facet_wrap(~ model)

Xem residual của từng line

sim3 <- sim3 %>% 
  gather_residuals(mod1, mod2)

ggplot(sim3, aes(x1, resid, colour = x2)) + 
  geom_point() + 
  facet_grid(model ~ x2)

Interactions (2 continuous)

mod1 <- lm(y ~ x1 + x2, data = sim4)
mod2 <- lm(y ~ x1 * x2, data = sim4)

grid <- sim4 %>% 
  data_grid(
    x1 = seq_range(x1, 5), 
    x2 = seq_range(x2, 5) 
  ) %>% 
  gather_predictions(mod1, mod2)
grid

## # A tibble: 50 × 4
##    model    x1    x2   pred
##    <chr> <dbl> <dbl>  <dbl>
##  1 mod1   -1    -1    0.996
##  2 mod1   -1    -0.5 -0.395
##  3 mod1   -1     0   -1.79 
##  4 mod1   -1     0.5 -3.18 
##  5 mod1   -1     1   -4.57 
##  6 mod1   -0.5  -1    1.91 
##  7 mod1   -0.5  -0.5  0.516
##  8 mod1   -0.5   0   -0.875
##  9 mod1   -0.5   0.5 -2.27 
## 10 mod1   -0.5   1   -3.66 
## # … with 40 more rows

Sử dụng seq_range() thay thế cho data_grid(). Thay thế cho việc sử dụng giá trị unique của x

Sử dụng các argument như

• Ví dụ

  seq_range(c(0.0123, 0.923423), n = 5)

  ## [1] 0.0123000 0.2400808 0.4678615 0.6956423 0.9234230

• pretty = TRUE

  seq_range(c(0.0123, 0.923423), n = 5, pretty = TRUE)

  ## [1] 0.0 0.2 0.4 0.6 0.8 1.0

• trim = 0.1, cắt 10% giá trị đuôi

  x1 <- rcauchy(100)
  print(x1)

  ##   [1]    7.17455599    0.51857630    1.50053731 -173.50882444   10.55462922
  ##   [6]   -0.12366276    0.03673073   -2.05993929    3.14016643    2.87335478
  ##  [11]    1.17228076   -1.78004277   -1.64912075    0.37686329   -0.93382796
  ##  [16]   -1.41491860   -1.27807239    0.60098207    1.15952963   -1.37894904
  ##  [21]    0.34351244  -15.64513911   -5.17460388   -0.67602844    4.66956926
  ##  [26]   -0.06596813   -0.42246329    0.66521295    0.45269676    0.46496412
  ##  [31]    4.80513029   -2.58605494   -5.32788465   -3.30938897   -1.68157543
  ##  [36]    3.77587123    0.58505480    0.30220331    0.43743750   -0.01776576
  ##  [41]   -0.09188035   10.58232172   -2.87903719   -0.84014695    8.91803085
  ##  [46]   -2.32556547   -4.61932400   -0.25750880   -0.31349105    0.29792883
  ##  [51]    2.83804175    2.08793754   15.70480489    1.30658619   -3.13945941
  ##  [56]    2.97463983   -1.68151356   -1.16133722   -0.71174224   -0.71249720
  ##  [61]    0.41676315    0.20540708    0.50853307   -0.62810086   -1.16936187
  ##  [66]    1.88497384    2.08855329   -0.96490480    0.15635953    2.08706652
  ##  [71]    1.62181829   -0.34904277    1.48756045   -0.19601606   -0.56082780
  ##  [76]   -0.18970034   -0.93878887  -25.24477495    8.03804283    5.57527568
  ##  [81]   -5.46078151   -0.41090739   -0.12711524    2.37673550    8.02653473
  ##  [86]   -3.11376245   -1.01051437    4.87194781   -1.47955698    0.41531512
  ##  [91]    1.29209637   -1.90283696   -6.13795523    0.25426971   -6.46131391
  ##  [96]   -1.28506105   -0.81624626   -0.03923384    0.64624844    2.45024668

  seq_range(x1, n = 5, trim = 0.1)

  ## [1] -5.494640 -2.114203  1.266235  4.646673  8.027110

vẽ

ggplot(grid, aes(x1, x2)) + 
  geom_tile(aes(fill = pred)) + 
  facet_wrap(~ model)

Không thấy rõ model nào fit tốt hơn

ggplot(grid, aes(x1, pred, colour = x2, group = x2)) + 
  geom_line() +
  facet_wrap(~ model)

ggplot(grid, aes(x2, pred, colour = x1, group = x1)) + 
  geom_line() +
  facet_wrap(~ model)

23.4.4 Transformations

log(y) ~ sqrt(x1) + x2 biến đổi thành log(y) = a_1 + a_2 * sqrt(x1) + a_3 * x2

Nếu transformation sử dụng dấu +, *, ^ hay -, ta phải bỏ vào I( )

y ~ x + I(x ^ 2) biến đổi thành y = a_1 + a_2 * x + a_3 * x^2

R sẽ tự động bỏ biến dư nếu viết dạng x + x

y ~ x ^ 2 + x sẽ thành y = a_1 + a_2 * x

df <- tribble(
  ~y, ~x,
   1,  1,
   2,  2, 
   3,  3
)
model_matrix(df, y ~ x^2 + x)

## # A tibble: 3 × 2
##   `(Intercept)`     x
##           <dbl> <dbl>
## 1             1     1
## 2             1     2
## 3             1     3

model_matrix(df, y ~ I(x^2) + x)

## # A tibble: 3 × 3
##   `(Intercept)` `I(x^2)`     x
##           <dbl>    <dbl> <dbl>
## 1             1        1     1
## 2             1        4     2
## 3             1        9     3

Chọn polynomial để smooth tốt

model_matrix(df, y ~ poly(x, 2))

## # A tibble: 3 × 3
##   `(Intercept)` `poly(x, 2)1` `poly(x, 2)2`
##           <dbl>         <dbl>         <dbl>
## 1             1     -7.07e- 1         0.408
## 2             1     -9.07e-17        -0.816
## 3             1      7.07e- 1         0.408

Khi sử dụng poly, ngoài range của data, polynomial sẽ tự động dự đoán positive and negative infinity. Có sử dụng natural spline

library(splines)
model_matrix(df, y ~ ns(x, 2))

## # A tibble: 3 × 3
##   `(Intercept)` `ns(x, 2)1` `ns(x, 2)2`
##           <dbl>       <dbl>       <dbl>
## 1             1       0           0    
## 2             1       0.566      -0.211
## 3             1       0.344       0.771

sim5 <- tibble(
  x = seq(0, 3.5 * pi, length = 50),
  y = 4 * sin(x) + rnorm(length(x))
)

ggplot(sim5, aes(x, y)) +
  geom_point()

Fit thử 5 loại spline model

mod1 <- lm(y ~ ns(x, 1), data = sim5)
mod2 <- lm(y ~ ns(x, 2), data = sim5)
mod3 <- lm(y ~ ns(x, 3), data = sim5)
mod4 <- lm(y ~ ns(x, 4), data = sim5)
mod5 <- lm(y ~ ns(x, 5), data = sim5)

grid <- sim5 %>% 
  data_grid(x = seq_range(x, n = 50, expand = 0.1)) %>% 
  gather_predictions(mod1, mod2, mod3, mod4, mod5, .pred = "y")

ggplot(sim5, aes(x, y)) + 
  geom_point() +
  geom_line(data = grid, colour = "red") +
  facet_wrap(~ model)

23.5 Missing value

df <- tribble(
  ~x, ~y,
  1, 2.2,
  2, NA,
  3, 3.5,
  4, 8.3,
  NA, 10
)

mod <- lm(y ~ x, data = df)

## Warning: Dropping 2 rows with missing values

mod <- lm(y ~ x, data = df, na.action = na.exclude)

nobs(mod)

## [1] 3